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Abstract 

In An asymptotic result on compressed sensing matrices, a new construction for com¬ 
pressed sensing matrices using combinatorial design theory was introduced. In this paper, 
we use deterministic and probabilistic methods to analyse the performance of matrices 
obtained from this construction. We provide new theoretical results and detailed sim¬ 
ulations. These simulations indicate that the construction is competitive with Gaussian 
random matrices, and that recovery is tolerant to noise. A new recovery algorithm tailored 
to the construction is also given. 
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1 Overview 


In 2006 Donoho, Candes, Romberg and Tao lamno] laid the foundations for a revolutionary 
new signal sampling paradigm, which is now called compressed sensing. Their key insight was 
that many real-world signals have the special property of being sparse - they can be stored 
much more concisely than a random signal. Instead of sampling the whole signal and then 
applying data compression algorithms, they showed that sampling and compression of sparse 
signals can be achieved simultaneously. This process requires dramatically fewer measure¬ 
ments than the number dictated by traditional thinking, but requires complex measurements 
that are ‘incoherent’ with respect to the signal. 

Candes, Romberg and Tao established fundamental constraints for sparse recovery; one 
cannot hope to recover A:-sparse signals of length N in less than 0{k\ogN) measurements 
under any circumstances. They then established that several classes of random matrices 
meet this bound asymptotically. Two important examples are the random Gaussian ensem¬ 
ble, which has entries drawn from a standard normal distribution, and the random Fourier 
ensemble, which consists of a random selection of rows from the discrete Fourier transform 

matrix mm- 

In m, a construction for compressed sensing matrices based on block designs and complex 
Hadamard matrices was introduced (see Construction [T] below). Here we add to the analysis 
of these matrices, both establishing new results on their compressed sensing performance, and 
providing details of extensive simulations. Our main results are the following. 

1. Theorem 10 of [3j establishes that an n x N matrix created via Construction [1] has the 
(Ai, t)-recovery property for all values of t < ^, where n is the number of rows in the 
matrix. In Section [3] we show that there exist vectors of sparsity at most \/^ which 
cannot be recovered. Recall that <1> has the (I"!, t)-recovery property if every t-sparse 
vector m is uniquely recoverable from its image <I>m by .^i-minimisation. 

2. In Section 0] we give a non-rigorous analysis of Construction [1] which suggests that 
sufficiently large matrices created via the construction allow recovery of most vectors of 
sparsity 0{-\/n logn). 

3. In Section[5]we provide detailed simulations which suggest that the recovery performance 
of matrices obtained from Construction [1] is comparable to that of the Gaussian ensemble 
for matrices with hundreds of rows and thousands of columns. Simulations also suggest 
that signal recovery is robust against uniform and burst noise. Algorithm running times 
are better than for Gaussian matrices. 

4. In Section [6] we propose a new algorithm for signal recovery, tailored to Construction 
m We show that for most vectors of sparsity at most y/n, this algorithm runs in time 
O(nlogn). 

2 Preliminaries 

For our purposes here, a Hadamard matrix of order r is an rxr matrix H such that each entry 
of FI is a complex number with magnitude 1 and HH* = (where H* is the conjugate 
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transpose). The prototypical example of such a matrix is the character table of an abelian 
group; for cyclic groups this matrix is often referred to as a discrete Fourier transform matrix 
or simply a Fourier matrix. When we restrict the entries to real numbers, we state explicitly 
that the matrix is a real Hadamard matrix. For background on Hadamard matrices, we refer 
the reader to |14j . 

A pairwise balanced design (V, B) consists of a set V of points and a collection B of subsets 
of V, called blocks, such that each pair of points is contained in exactly A blocks for some 
fixed positive integer A. If n = |F| and K is a finite set of integers such that \B\ G K for each 
B ^ B, we use the notation PBD(u, IT, A). We denote the maximum element of K by 
and the minimum element of K by A"min- For each point x GV , the replication number r^ 
of X is defined by r^ = \{B G B : x G B}\. A set of points in a PBD is an arc if at most two 
of the points occur together in any block. 

A PBD(n, {/c}. A) is a balanced incomplete block design, denoted by BIBD(n,/c. A). Obvi¬ 
ously, all points of a BIBD must have the same replication number. This paper is devoted to 
the study of matrices created via Construction [H which generalises a construction from m- 

Construction 1 ([!]). Let (V,B) be a PBD(n,Ar, 1) with n = \B\ and N = Ylxev’’’^- Then 
d* is the n X N matrix constructed as follows. 


• Let A be the transpose of the incidence matrix of {V,B)-, rows of A are indexed by 
blocks, columns of A by points, and the entry in row B and column x is 1 ii x G B 
and 0 otherwise. 


• For each x gV ,\ei Hx be a (possibly complex) Hadamard matrix of order . 


For each x gV , column x of A determines columns of <I>; each zero in column x is 
replaced with the 1 x row vector (0,0,... , 0), and each 1 in column x is replaced 
with a distinct row of —^Hx- 

y/Fx 


Remark 2. Suppose that an nxN matrix is created via Construction [T] using a PBD(u, K, 1) 
V with replication numbers ri,...,and block sizes ki,... ,kn. We briefly discuss the 
behaviour of n and N in terms of the parameters of P in the case where v is large and 
Kmin ~ Aimax, that is, there exists a constant a not depending on v such that iLmax < aiLmin ■ 
Let r = + • • • + r^) and k = ^{ki + ■ ■ ■ + kn) be the average replication number and 

average block size respectively. Obviously ATmin, Lfmax ~ k, and standard counting arguments 
for block designs (see Chapter II of [T] for example) yield 


N = rv = kn 

( 1 ) 

kr ^ V 

( 2 ) 

1 2 2 
k n ^ V 

( 3 ) 


Combining ([2]) and ([3]), we see that n ~ r^. It is known (see [9]) that any non-trivial PBD 
has at least as many blocks as points. So n > v, and hence k grows no faster than 0(\/v) by 
©. Two extremes can be distinguished. When k is constant, n ~ r ~ -^/n and N ^ r'^ ^ n 
- this is the case considered in [3j. When k V^, y r'^ n and N ^ r^ ^ (this occurs 
when P is a projective plane, for example). It is easy to see that the growth rates of v and 
N are bonnded by these two extremes. 
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Except when n G {1,2}, there do not exist real Hadamard matrices of order n for n ^ 0 
(mod 4). If complex Hadamard matrices are employed in Construction [H the resulting com¬ 
pressed sensing matrix also has complex entries. If real matrices are required for a particular 
application (for example, linear programming solvers generally require the ordered property of 
the real numbers), the isometry in Lemma [3] can be exploited. While this result is elementary 
and well known, we first became aware of its use in this context in [3]. 


Lemma 3. The map 


a + ih\-^ 


a b \ 
—b a ) 


is an isometry from C to M 2 (M). It extends to an isometry from Mn,N{C) to M2n,2N0^)- 


This technique allows a complex compressed sensing matrix to be converted into a real 
one, at the expense of doubling its dimensions. If a complex matrix recovers t-sparse vectors, 
so too does the corresponding real matrix. (But while a t-sparse complex vector is in general 
2t-sparse, we cannot, and do not claim to, recover arbitrary 2t-sparse real vectors.) Where 
confusion may arise, we specify the field over which our matrix is defined as a subscript; thus 
if ‘he is n X N, <I>]k is 2n x 2N. 

We require a result concerning the sparsity of linear combinations of rows of a complex 
Hadamard matrix, which is essentially an uncertainty principle bounding how well the stan¬ 
dard normal basis can be approximated by the basis given by the (scaled) columns of a 
Hadamard matrix. In particular, if m G C” admits an expression as a linear combination of 
u columns of a Hadamard matrix, and as a linear combination of d standard basis vectors, 
then du>n. 


Lemma 4 ([18], Lemma 2, cf. |15] . Lemma 14.6). Let H be a complex Hadamard matrix of 
order n. If m is a non-zero linear combination of at most u columns of H, then m has at 
least [-] non-zero entries. 

The bound in LemmaHjis sharp. Let LI be a Fourier matrix of order u, and suppose that 
u divides v. Then there exist u columns of H containing only roots of unity and the 
sum of these columns vanishes on all but - coordinates. At the other extreme, a non-trivial 
linear combination of any u of the columns of a Fourier matrix of prime order vanishes on at 
most u — 1 coordinates. 


3 Upper bounds on performance 

The spark of a matrix <h is the smallest non-zero value of s such that there exists a vector 
m of sparsity s in the nullspace of <!>. In this section we provide upper and lower bounds on 
the spark of a matrix obtained from Construction [TJ 

The following well-known result bounds the (Ii, t)-recoverability of a matrix in terms of 
its spark. 

Proposition 5. If the spark of a matrix ^ is s and <I> has {ii,t)-i"Gcoverability, then t < §. 
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Proof. Let be a matrix with spark s and let m be an element of sparsity s in the nullspace 
of Write m = mi + m 2 where mi has sparsity [|J and m 2 has sparsity l"^]. Then 


<hmi + ‘hm2 = 0. 

So ‘h(m 2 ) = <h(—mi), and one of —mi or m 2 is not recoverable. Thus ‘h does not have 
(£ 1 , [|])-recoverability and the result follows. □ 

We can now provide upper and lower bounds on the spark of <h. 

Proposition 6. Let P be a PBD(u,iL, 1) whose smallest replication number is ri and let 
be a matrix obtained from Construction [7] using P. Then the spark s of ^ satisfies ri < s. 

Proof. Suppose that m is in the nullspace of <h. For each i G u}, let m* be 

the vector which is equal to m on those coordinates corresponding to point i of P and 
has each other coordinate equal to 0. Let T = {i : supp(mj) / 0 }, and t = \T\. 
Now, $mj = ^>m = 0 and, because any two points of P occur together in ex¬ 

actly one block, |supp(<hmj) n supp($mj)| < 1 for any distinct i,j € T. Thus it 
must be that |supp(^>mi)| < t — 1 for each i G T. By Lemma 01 this implies that 
|supp(mi)| > for each i ^ T. So, because ri > ri for each i G T, we have that 
|supp(m)| = JfieT |supp(mi)l > ^ > ri. □ 

Proposition 7. Let P be a PBD(u,iL, 1) whose two smallest replication numbers are ri and 
^2 (possibly ri = r^) and let ^ be a matrix obtained from Construction{l\ using P. Then the 
spark s of ^ satisfies s < ri + r 2 . 

Proof. Let pi and p 2 be points of P with replication numbers ri and r 2 respectively, and 
let <hi and ^2 be the submatrices of ‘h consisting of the columns corresponding to pi and p 2 
respectively. We establish the result by showing that the columns of the block matrix [<hi <^ 2 ] 
are linearly dependent. 

There is a unique row p of the matrix [^>i <^ 2 ] that contains no zero entries (corresponding 
to the unique block that contains pi and P 2 ). For i G {1,2}, let be the matrix obtained 
by scaling the columns of so that all entries in p have the value ^ if i = 1 and — ^ if 
i = 2. For i G {1,2}, because the rows of <l>i are pairwise orthogonal, the rows of are too. 
In particular, row t of is orthogonal to each other non-zero row of and thus the sum 
of any row of other than row t is zero. Clearly, the sum of row p of <hi is 1 and the sum 
of row /9 of <^2 is — 1. Thus the columns of add to the zero vector, completing the 

proof. □ 

The bound of Proposition [7] is sharp. Let P be a PBD(u,iL, 1) with every replication 
number prime, and let be a matrix obtained from Construction [T] taking {V,B) to be P 
and Hx to be a Fourier matrix of order for each x gV . Using the fact that a non-trivial 
linear combination of u columns of a prime order Fourier matrix vanishes on at most u — 1 
coordinates, it can be shown that the spark of $ is exactly ri -|- r 2 . We say that a Hadamard 
matrix is optimal if no linear combination of k rows has more than k zero entries. By this 
criterion the Fourier matrices of prime order are optimal. The next remark shows that, in 
general, the spark of depends both on the choices of the design and the Hadamard matrices 
used in Construction [TJ 
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Remark 8. Let P be a PBD(?;,iL, 1) with every replication number equal to some r = 0 
(mod 4), and let be a matrix obtained from Construction [1] taking (V,13) to be P and 
Hx to be a real Hadamard matrix of order r for each x G V. Denote by 1 a vector of Is of 
length I. Consider the submatrix <!>' of 4> consisting of the columns corresponding to three 
points of an arc of P. Then by scaling and reordering columns of we can obtain a matrix 
that has a submatrix of the form 




The vector m = (0,1,0,-1,0,1)T is in the nullspace of This is easily verified for 
the displayed rows, and follows for the remaining rows from the orthogonality of Hadamard 
matrices (the vector m is a linear combination of the displayed rows of the Hadamard matrix 
when restricted to any set of columns corresponding to a point). It follows that the nullspace 
of <I> contains elements of sparsity ^. 

Real Hadamard matrices are never optimal, by Remark [51 The obvious questions here 
concern the existence of optimal or near optimal complex Hadamard matrices; some discussion 
is contained in [18] . 

Now we combine our results so far to determine an upper bound for the {ii , t) -recoverability 
of 4>. Recall that MIP-based results are sufficient conditions for (^i, t)-recoverability, but are 
not necessary in general. Thus, while the Welch bound forms a fundamental obstacle to using 
MIP to prove that vectors of sparsity exceeding ^ sampled by an n x iV matrix can be 
recovered, it does not guarantee that such vectors cannot be recovered. We show that there 
exist vectors of sparsity at most which cannot be recovered by matrices obtained via 

Construction [1] 


Proposition 9. Let V he a PBD(u,iL, 1) such that Rlmax < — 1) cind the two 

smallest replication numbers of P are ri and r 2 (possibly ri = r 2 ). Let <I> be an n x N 
matrix obtained from Construction{l\ using T>. If ^ has {ii,t)-recoverability, then t < y/^. 


Proof. By Proposition [71 the spark s of 4> is at most ri -|-r 2 . The replication number of any 


point is at most , so 

-^min J- 


s < 


2iv - 1) 

K — 

■tkmin 


Because n is the number of blocks in P and R'max < V^{Kra\n — 1), we have 


n > 


v{v — 1) 


7ilmax(7iimax 1) 


{V - 1)2 

> - — > 

- 


{v - 1)^ 


> 


2{K^^ - 1)2 - 8 


and thus s < 2y/^. The result now follows by Proposition |5l 


□ 


Combining Proposition [9] with Theorem 10 of [1] we can specify the (^i, t)-recoverability 
of a matrix obtained from Construction [1] to within a multiplicative factor of 4\/2. 

Theorem 10. Let T> be a PBD(u,iL, 1) such that Rlmax < ■\/2(7Cmin~l) cind the two smallest 
replication numbers of P are ri and r 2 (possibly ri = r 2 )■ Let <I> be an n x N matrix 
obtained from Construction [i| using P and let t* be the greatest integer t such that <I> has 
{ii,t)-recoverability. Then t* = c^fn for some c in the real interval [l,-\/2]- 


6 









A matter of practical concern is whether, despite the presence of vectors of a certain 
sparsity that cannot be recovered, one can nonetheless recover many vectors of larger sparsity. 
We consider a sample computation. 

Example 11. An arc of maximal size in PG(2,(/) is called an oval. When q is odd, a 
celebrated result of Segre shows that an oval is the set of points on a conic and has size q + 1 
(see Chapter 8 of m)- A convenient construction for ovals is as the negation of a Singer set 
(Proposition VII.5.12 of [T] ). The projective plane PG(2,7) corresponds to a BIBD(57,8,1). 
Remove eight points of an oval to produce a PBD(49, {6, 7,8}, 1). Apply Construction [1] with 
(VjB) as the PBD using Fourier matrices of order 8, to create a 57 x 392 matrix <I>c. Apply 
Lemma[3]to create a 114 x 784 matrix <h]R. The spark of ‘he is at most 16 by Proposition [71 
So, by Proposition ISl there exist 8-sparse vectors that cannot be recovered by <hc ~ this is in 
agreement with the bounds 1 < t < 10 obtained from Theorem [TOl (with n = 57). 

Figured] shows the number of successful recoveries of vectors of sparsity t sampled by <I>ir 
out of 1,000 attempted recoveries using an OMP algorithm [T7] for each t € {1,2,... ,50}. 
(See Section [5] for details of all simulations.) 



Figure 1: Signal recovery for a 114 x 784 matrix obtained from Construction [T] 


Example dH suggests that although there are vectors of sparsity less than ^/^ which 
cannot be recovered by an n x V matrix <I> obtained from Construction dl such vectors are 
rare. In fact, the typical performance of <h is much better. In the next section we explore 
heuristic arguments for why this should be so. 

4 Heuristic arguments for small sparsities 

Throughout this section we suppose that <I> is an n x V matrix obtained from Construction 
dl using a PBD D with all replication numbers equal to r. As in Example [TTl simulations 
suggest that although there are vectors of sparsity 2r in the nullspace of <I>, and hence there 
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exist vectors of sparsity r which cannot be recovered, such vectors are rare. In this section, 
we give a heuristic justification, inspired by techniques from random matrix theory. The 
exposition in [20] is relevant here. 

4.1 Random matrix theory 

Let M be an n X n Hermitian matrix. Then the eigenvalues Ai(M) > • • • > A„(M) of M 
are real and the function Lm{o) = ^\{i ■ \{M) < a^/n}\ defines the density function of a 
(discrete) probability distribution on M. 

One of the central problems of random matrix theory is to describe the function Lm when 
M is drawn from some specific class of matrices. The fundamental result, due to Wigner, 
concerns real symmetric nxn matrices with all lower triangular entries drawn independently 
from a Gaussian (0,1)-distribution. His main result was that as n ^ oo, Lm converges to 
the semi-circle distribution: 

L{s) = ^j Vl - ds. 

Wigner originally demonstrated pointwise convergence and later contributors obtained 
stronger results, weakening assumptions on the probability distribution of the entries of M 
and establishing convergence of measure. We require only one of the bounds on the eigenvalues 
of M, though there is ample scope for further application of random matrix theory to the 
analysis of deterministic matrix constructions. The operator norm of a Hermitian matrix M , 
denoted ||M||op, is the maximum absolute value of an eigenvalue of M. Following |20] we 
will say that an event E depending on a parameter t occurs occurs with high probability if 
¥{E) > 1 —o(l) and with overwhelming probability if, for every fixed H > 0, P(£') > l — CAt~^ 
for some Ca not depending on t. 

Lemma 12 1[20]. Corollary 2.3.6). Let M = (mij) be a random 2tx2t Hermitian matrix such 
that the entries for 1 < i < j < 2t are jointly independent with mean 0 and magnitude 
uniformly bounded by 1. Then there exist absolute constants C,c > 0 such that for any real 
number A> C, 

P ^||M||op > < Cexp{—cAt). 

In particular, for large t, ||M||op = 0{^/i) with overwhelming probability. 

4.2 RIP bounds and signal recovery 

We recall that <h has the restricted isometry property with constants t, 6, abbreviated RIP(t, <5) 
if and only if for each t-sparse vector m, 

(1 - d)\\m\\l < \\^m\\l < (1 + S)\\m\\l. 

This is equivalent to the statement that, for each set S of at most t columns of <h, the 
eigenvalues of the matrix all lie in the interval [1 — <5,1 + 5], where is the submatrix 

of 4> containing only the columns in S. 

RIP conditions have been used extensively to provide sufficient conditions for {li,t)- 
recoverability [5l|6|. For our purposes it is enough to note that if satisfies the RIP(2t, ^/2—l ), 
then $ has (t'l, t)-recoverability. In the next subsection, we develop a simple model for signal 
recovery, which suggests that <I> recovers vectors of sparsity 0{^fnlogn) with high probability. 




4.3 A heuristic model 


We begin by developing a model for signals m with the property that all non-zero coordinates 
occur in columns of corresponding to different points of the design D. Recall that we are 
assuming that all points in D have equal replication number r. 

Let S be a set of 2t columns of <h, no pair corresponding to the same point, and denote 
by ‘hs' the submatrix of <1> consisting of the columns in S. With these assumptions, is 

of the form I + where is a Hermitian matrix with zero diagonal, and all off-diagonal 
entries of magnitude ^ . We lack information on the phase of the entries of 4 ' 5 . Consider the 
space of matrices 4^5 as S' varies over all sets of 2t columns of <1> with no pair corresponding 
to the same point. Our heuristic is that a typical element of this space behaves as the matrix 
^M, where M is a 2t x 2t random Hermitian matrix with zero diagonal in which the strictly 
upper triangular entries have magnitude 1 and uniformly random phase. 

Now, observe that the eigenvectors of / -|- 'Lg are those of 4 ' 5 . By Lemma [T^ for large t, 
||.^L||op = 0{\/t) with overwhelming probability. Thus, provided that t = o(r^) as r —>■ oo, 
the eigenvalues of = I+ s are all arbitrarily close to 1 with overwhelming probability. 

In particular, our heuristic suggests that 4' recovers any vector with the property that its 
support intersects each point of the design in at most one column with high probability. 

Allowing d non-zero entries in S to be labelled by the same point of V introduces 2 ( 2 ) 
off-diagonal zero entries in T 5 . The eigenvalues of a matrix are continuous functions of the 
matrix entries (via the characteristic polynomial). So provided that the total number of zeroes 
introduced is not excessive, the analysis of the heuristic continues to hold. Thus, it seems 
reasonable to suppose that our heuristic is valid for signals of sparsity at most r log r with at 
most log r non-zero coordinates corresponding to any point of the design. We call such signal 
vectors suitable. 

What we have obtained so far is not a uniform recovery guarantee, however; it does 
not follow that <I> recovers all suitable vectors with high probability. We need a slightly 
more careful analysis to obtain a uniform recovery guarantee. Taking t = rlogr and 
A = ^(2 — \/2)r^/^(logr)“^/^ in Lemma fT^ we obtain 

P (i||M||op > V 2 - 1 ) < 

for absolute constants C and c. Now, the total number of subsets of size 2rlogr of a set of 
size N = vr = O(r^) is bounded above by . So, because r^/^(log 

grows faster than 6 r(logr)^, we have that P(i||M||op > \/2 — 1) = . It follows 

by the union bound that for sufficiently large r the eigenvalues of all lie in the interval 
[2 — ^/2, \/2] with high probability. (Recall that r is the replication number of the design used 
in Construction [ 1 ] to obtain 4> and that the order of M is a function of r.) 

Consequently, the heuristic suggests that there is a high probability that, for large r, 4> 
satisfies the RIP(2r logr, \/2 — 1) on all suitable sets of columns, and so all suitable vectors 
are recoverable. By Proposition [9] and Lemma [5] we know that there exist vectors of sparsity 
r that are not recoverable, but these are far from suitable. 

This analysis suggests that for large r, <I> recovers all suitable vectors of sparsity yjn log n , 
with high probability. This result is in excess of the square-root bound. The probability that 
4> fails to recover any particular random signal vector is exponentially small, being essentially 
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the probability that many non-zero coordinates of the signal are concentrated in columns 
labelled by a small number of points. 

5 Simulations 

In this section we describe the results of extensive simulations. All simulations are performed 
with real-valued matrices and signal vectors. If our construction yields a complex matrix, we 
apply Lemma [3] to obtain a real one. To illustrate the flexibility of Construction [1] we derive 
input PBDs from BIBDs in various ways. 

Our matrices are stored as an array of real floating point numbers of some fixed accuracy 
1. The entries typically consist of numbers of the form cos(2s7r/r) or sin(2s7r/r) where r is 
a replication number of the design used to create the matrix and s G {0,... ,r — 1}. Our 
simulations are performed as follows. 

• Vectors of sparsity t are constructed by choosing t coordinate positions uniformly at 
random, and populating these locations either with a value from a uniform distribution 
on (0,1) (in matlab) or a uniform distribution on : 1 < * < 100} (in MAGMA) and 
then normalising. In all discussion of simulations, when we refer to a t-sparse vector we 
mean that it has exactly t non-zero coordinates. 

• We invoke the standard implementation of one of the recovery algorithms from Section 
[Q to obtain a vector m such that ~ $m. Only and are provided to the 
solver; the sparsity of m is not supplied. 

• The other parameters in our test are the precision to which real numbers are stored 
(typically 25 decimal places) and the recovery error allowed, e. A trial is counted a 
success if |m — m| < e, and a failure otherwise. We report the proportion of successes 
at sparsity t as a proxy for compressed sensing performance. 

Some general comments apply to all of our simulations. In essence choosing a design 
and Hadamard matrices for the construction is a multi-dimensional optimisation, probably 
application specific. After fixing a matrix ^>, one chooses the precision, the maximal entry 
size in a random vector and the error allowed in recovery. These all impact performance, both 
proportion of correct recoveries and recovery time. 

The recovery process is not very sensitive to the allowable recovery error, in the sense that 
if the algorithm converges to the correct solution it generally does so to machine tolerance. 
Taking sparse vectors to be binary valued greatly improves recovery, but we focus on the 
model given. Numerical instability results if the precision to which real numbers are stored is 
too small; this is determined experimentally. 

5.1 Recovery algorithms 

We compare the performance of two different well-established recovery methods when applied 
to signal vectors that have been sampled using matrices created via Construction [TJ 

Firstly, we employ naive Linear Programming (LP); we consider implementations in both 
MAGMA and matlab, mm- We include implementations in two different systems to demon¬ 
strate the potential differences between solvers, which can be significant. 
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Secondly, we employ matching pursuit, which is a greedy algorithm for producing sparse 
approximations of high dimensional data in an over-complete dictionary. A well-known im¬ 
plementation for compressed sensing is CoSaMP m- We use the implementation of the basic 
OMP algorithm from m rather than the more specialised CoSaMP algorithm. It is well 
known that the worst case complexity of the simplex method in linear programming is expo¬ 
nential in the number of variables. However the expected complexity for typical examples is 
0{N^). In comparison, the complexity of CoSaMP is 0(iV log^A^). 

We created a 57x456 matrix <hc using ConstructionUl taking (P, B) to be a BIBD(57,8,1) 
corresponding to a PG(2,7) and Hi,..., Hy to be Fourier matrices of order 8. We then 
applied Lemma [3] to obtain a 114 x 912 matrix <I']r. The allowable recovery error was 10“® 
for all algorithms. In Figure[2l we record the number of successful recoveries for 1000 randomly 
generated vectors sampled by ‘I'k for each algorithm and each sparsity. 

Applying Theorem [TO] shows that <hc allows recovery of all vectors of sparsity at most 
= 2, and provides no stronger guarantee. Furthermore, <l>c contains real rows so, 
as in Remark [8l there exist 12-sparse vectors in the nullspace of 4>c, and hence there exist 
6-sparse vectors whose recovery is not permitted by 4>c ■ So the recovery performance of <h]K 
is quite striking. 



Figure 2: Comparison of recovery algorithms 

While there are clear differences in performance among these algorithms, they appear 
to behave in a broadly similar fashion. In general, running times for OMP are an order 
of magnitude faster than the matlab linear programming solver. The MAGMA solver has 
intermediate runtime. 

5.2 Comparison with Gaussian matrices 

Gaussian matrices are the de facto standard against which other matrix constructions 
are measured in compressed sensing. The projective plane PG(2,11) corresponds to a 
BIBD(133,12,1). Removing two blocks and all points incident with either block produces 
a PBD(110, {10,11}, 1) in which all points have replication number 12. We applied Con- 
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structiondl taking {V, B) to be this PBD and Hi ,..., to be Fourier matrices of order 12, 
to obtain a 131 x 1320 matrix <I>c. We then applied Lemma[3]to obtain a 262 x 2640 matrix 

We compare the compressed sensing recovery performance of to that of a Gaussian 
ensemble with the same numbers of rows and columns (using the OMP algorithm). Both linear 
programming and OMP ran faster by an order of magnitude for 4 >]r than for the Gaussian 
ensemble. Our results are given in Figure [3l 



Figure 3: Gomparison with Gaussian ensemble 


5.3 Factors affecting algorithm performance 

In this section we discuss a number of factors that influence recovery performance. Specifically 
we consider the presence of noise in received signals, signed signal vectors (until now our signal 
vectors have had positive coordinates), and the effect of using different Hadamard matrices 
in Gonstruction [TJ In all cases, we find that the construction is robust and reductions in 
performance are not substantial. 

5.3.1 Signed signal vectors 

The linear programming solver in MAGMA requires variables to be positive (or at least 
greater than some bound). We use a standard trick to allow negative entries in x. For each 
variable Xi in the original problem, we introduce a pair of variables, xf and x^. Then we 
replace each appearance of Xi with xf — x~. 

When a solution has been found, we interpret xf as a positive number, and x^ as a 
negative one. This has the disadvantage of doubling the number of variables in the linear 
program. This results in longer run times and slightly lower performance, but the results of 
the simulations are broadly comparable to those in the positive case. 
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5.3.2 Noise 


We consider both uniform and burst noise. In each case we consider positive and signed noise 
vectors. In general, recovery in the presence of noise is robust. Testing all algorithms in 
all regimes would produce an overabundance of data, so we give only a single representative 
example of our simulations. 

The projective plane PG(2,11) corresponds to a BIBD(133,12,1). Remove twelve points 
of an oval to produce a PBD(121, {10,11,12}, 1). Apply Construction [T] taking Hi ,..., 
to be Fourier matrices of order 12, and then apply Lemma [3] to construct a compressed 
sensing matrix <1 >r with 266 rows and 2904 columns. To examine the performance of <1 >]r 
in the presence of uniform positive noise, we construct noise vectors with entries uniformly 
distributed in (0,1) and scale the vector to some predetermined .^ 2 -norm. We then compare 
the signal vector m (of £2 -norm 1) to the solution returned by the matlab linear programming 
solver for 4>R(m + e). We count a recovery as a success if the reconstruction error is below 
10“®. Table [U summarises our results. 


Sparsity 

1 ' 2 -norm of noise vector 

0 

10"^^ 

10-iu 

10"^ 

2 X 10"^ 

30 

100 

99 

98 

79 

66 

35 

100 

100 

97 

79 

69 

40 

100 

100 

91 

77 

49 

45 

97 

93 

88 

62 

27 

50 

87 

79 

69 

33 

5 

55 

61 

56 

30 

14 

2 

60 

27 

22 

22 

5 

0 


Table 1; Number of successful recoveries out of 100 for signals of various sparsities and for 
different noise levels. 


Recovery decays gracefully in the presence of noise, particularly when the sparsity of the 
signal is not close to the limit of what can be recovered. 

5.3.3 Choice of Hadamard matrix 

We illustrate the effect of the choice of Hadamard matrix with a small example. We take a 
BIBD(25,3,1) obtained from the Bose construction, which has replication number r = 12 
(see p.25 of El)- The 100 X 300 matrices <l>c and are obtained via Construction [H 
taking {V,B) to be this BIBD and Hi,..., to be real Hadamard matrices of order 12 
and Fourier matrices of order 12, respectively. Lemma [3] is then applied to both matrices 
to obtain 200 x 600 matrices ‘I'k and (While is already real, this application of 

Lemma [3] allows a direct comparison between constructions. It has no effect on the sparsity 
of vectors recovered.) Recovery performance varies by less than 2% and runtime by less than 
6% with the OMP algorithm. Such variation could be caused simply by random fluctuations. 
On the other hand, the differences in performance in the MAGMA LP-implementation are 
substantial. We record them in Table El 
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Sparsity 

Real matrix 

Fourier matrix 

No. successes 

Avg. time 

No. successes 

Avg. time 

56 

100 

4.8 

100 

(crashed twice) 

58 

100 

4.5 

99 

14.5 

60 

98 

4.4 

96 

14.8 

62 

100 

4.6 

98 

15.1 

64 

93 

4.8 

99 

15.6 

66 

96 

4.8 

97 

15.9 

68 

94 

4.9 

96 

16.3 

70 

87 

5.0 

99 

16.8 


Table 2: Number of successful recoveries out of 100 and the average recovery time in seconds 
for real and Fourier matrices. 


With the MAGMA LP-implementation, it would appear that <l>{g produces better recovery 
at the cost of increased runtime and the risk of numerical instability. Obviously the choice of 
matrix and recovery algorithm would depend on the particular application. 

6 An efficient algorithm for sparse recovery 

In this section we describe and investigate a new algorithm for signal recovery, tailored specif¬ 
ically for an n X N compressed sensing matrix created via Construction [TJ It is designed to 
recover vectors of sparsity at most 0(^/n) and is not expected to be competitive with LP 
or CoSaMP at large sparsities. The algorithm exploits the structure of matrices created via 
Construction [1] in order to achieve efficiency in both running time and storage requirements. 
Under certain assumptions, it can be shown to run successfully in time 0{N\ogn) and space 
0{'n?). (Such thresholding based algorithms are one of the main algorithmic approaches to 
compressed sensing; see Chapter 3 of [E].) Throughout this section we employ the notation 
that we introduce in describing the algorithm. 

Algorithm 13. Suppose that <I> is a matrix created via Construction [U using a design (V,13) 
with replication numbers ri < • • • < and Hadamard matrices Hi ,..., . For each i G U, 

let <I>j be the submatrix of consisting of the columns corresponding to point i. We store 
only the following information. 

• For each i € U, a list of the rows of that are non-zero, and a record of which rows 
of Hi are located there. 

• A copy of each Hadamard matrix Hi. 

Suppose that $ is used to sample some N x 1 message vector m. Our algorithm runs as 
follows. 

1. Construct an initial estimate. For each i £ V, we take the set of rows in which <I>j is 
non-zero, take the rj x 1 vector yi of the corresponding received samples, and compute 
ifii = -^H*yi. Concatenating the vectors thi over all i G U we construct an initial 
estimate m for the signal vector m. 
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2. Guess the signal coordinates. For some prespecified |5| < ri, let S be the index set of 

the I S'! coordinates of m of greatest magnitude. For each i € let qi be the number 
of columns indexed by S that are in d>j, and let V' = {i : Qi > 1} ■ 

3. Identify nncontaminated samples. For each i £ V', we find a set Qi of qi rows of $ 

that are non-zero in but zero in for each j £ V' \ {i} (such rows correspond to 
blocks that contain i but no point in and, because < \S\ — qi < ri — qi, 

at least qi such blocks exist). 

4. Recover the signal. For each i £ V', we hnd the coordinates in the qi positions indexed 
by S that correspond to columns in by solving the qi x qi linear system induced by 
those columns and the rows in Qi. 

Remark 14. The record of the non-zero rows of <I> can be derived easily from the incidence 
matrix of {V,B). For many designs (for example, projective planes, designs with cyclic auto¬ 
morphisms, and so on) this information need not be stored explicitly. Even if the information 
is stored explicitly, the space used is O(n^). A similar observation holds for the Hadamard 
matrices; Fonrier or Paley matrices need not be stored explicitly. With appropriate design 
choices, the storage space required can be logarithmic in n. 

As long as S contains the positions of all the non-zero elements, the algorithm hnds the 
right solution. In our analysis of this algorithm we focus, for the sake of simplicity, on the 
case in which all the replication numbers of (y,B) are equal. We first show that our estimate 
for a signal coordinate differs from the actual value by an error term introduced by the other 
non-zero signal coordinates. From this it is easy to show that Algorithm 1131 recovers a signal 
provided that no non-zero signal coordinate has magnitude too small in comparison with the 
£i-norm of the signal. 

Lemma 15. Let ^ be a matrix created via Construction [I] using a PBD V all of whose 
replieation numbers are equal to some integer r. Suppose that is used to sample a signal 
m of sparsity at most r and that Alaorithm\13\ is applied to find an estimate rh for m. Let 
X be the estimate for a coordinate x of m that corresponds to a point i of V. Then 

(i) X = X + ^{h\X\ -|- • • • -b hdXd) where xi,... ,Xd are the non-zero coordinates of m that 
do not eorrespond to point i of V and hi,..., are complex numbers of magnitude 1; 
and 

(a) |x — x| < y ||m||i. 

Proof. Let xi,...,Xd be the non-zero coordinates of m that do not correspond to point i 
of V and let Wi be the vector (xi,... ,Xd)^ ■ Let be the r x I vector consisting of the 
coordinates of m corresponding to point i of P and let rhi be the estimate for mi. Note 
that yi = -^Himi -£ AiWi for some r x d matrix Ai such that each column of Ai contains 
exactly one complex number of magnitude 1 and every other entry of Ai is a 0. So 

rhi = lH*Himi + -^H*AiWi = mi + ■^{H*Ai)wi. 

Thus, because every entry of H*Ai is a complex number of magnitude I, 

x = x + \{hixi H-h hdXd) 
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where hi,... ,hd are complex numbers of magnitude 1. Further, 


x-x\ = \^{hixi H-h hdXd)\ < i(|xi| H-h \xd\) < 


□ 

Corollary 16. Let ^ be a matrix created via Construction^ using a PBD T> all of whose 
replication numbers are equal to some integer r. Suppose that is used to sample a signal 
m of sparsity t <r. If each non-zero coordinate of m has magnitude at least |||m||i, then 
Algorithm{T^with l^l > t recovers m. 

Proof. Let xq and xi be the estimates obtained by Algorithm [13] for two coordinates xq and 
xi of the signal such that xq = 0 and xi / 0. It suffices to show that |xo| < |xi|. By 
Lemma [T5lf hi. |xo| < ^||m||i. From our hypotheses |xi| > y||m||i and so by Lemma [T5jii) , 
|xi| > |xi| - i||m||i > y||m||i. Thus, |xo| < |xi|. □ 

If we assume that the non-zero signal coordinates have uniformly random phase, we can 
improve substantially on Corollary[T6| Also, if we assume that the support of the signal is cho¬ 
sen uniformly at random, we can establish that, for large n, the algorithm runs in 0{N \ogn) 
time with high probability. We employ a simple consequence of Hoeffding’s inequality. 

Lemma 17. If zi,... ,Zn are independent complex random variables with uniformly random 
phase and magnitude at most 1, then for each positive real number c, 

H-h 2 n| > c) < 4exp • 

Proof. If \zi + ■ ■ ■ + Zn\ > c, one of the real or imaginary parts of zi + ■ ■ ■ + Zn must have 
magnitude at least . Now apply Hoeffding’s inequality separately to the real and imaginary 
parts of zi + ■ ■ ■ + Zn, noting that the expected value is zero in each case. □ 

Theorem 18. Let d> be an n x N matrix created via Construction{l\ using a PBD D all of 
whose replication numbers are equal to some integer r. Suppose that d> is used to sample an 
r-sparse signal m. 

(i) If the non-zero components of m are independent random complex variables with uni¬ 
formly random phase and magnitude in the interval [r~^^^, 1] for a positive constant e, 
then for large n Alaorithm\lS\ with I^I = r recovers m with high probability. 

(ii) If the support of m is chosen uniformly at random, then for large n the algorithm runs 
in 0{Nlogn) time with high probability. 

Proof. Recall that n ^ r'^ and suppose that r is large. 

We hrst prove (i). For brevity, let a = Suppose that the non-zero components of 

m are independent random complex variables with uniformly random phase and magnitude 
in the interval [o, 1]. Let x be a coordinate of m and let x be our estimate for x. By Lemma 
m i), X — X = -^(hixi hdXd) where xi ,... ,Xd are the non-zero coordinates of m that 

do not correspond to point i of D and hi,...,hd are complex numbers of magnitude 1. 
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Because xi,... ,Xd are independent random complex variables with uniformly random phase 
and magnitude at most 1, so are hixi, ..., h^Xd- Thus, by Lemma [T71 

P + • • • + hdXd\ > f) < 4exp = 4exp 

Thus, using the facts that d < r and that grows faster than log N (by Remark [2l 
N = 0{r^) ), it can be seen that P(|x — x| > |) = o{^). So it follows from the union bound 
that with high probability it is the case that |x —x| < | for each coordinate x of m. Because 
|x| > a for each non-zero coordinate x of m, this implies that our estimates for non-zero 
coordinates of m have magnitude greater than our estimates for zero coordinates of m with 
high probability. Then Algorithm 1131 with l^l = r recovers m. 

We now prove (ii). Note that l^l = r, that rv = N and that r ~ ^/n (see Remark [2]). 
The first step of the algorithm is essentially performing a Hadamard transform of order r 
for each point, and so can be accomplished in 0(ur logr) = 0{N \ogn) time. The second 
step is essentially a sorting operation and can be accomplished in 0{N) time. The third step 
can be accomplished by first creating a list of blocks that intersect exactly one point in V' 
(by examining r blocks for each of the at most \S\ points in V'} and then partitioning this 
list into the sets Qi. Becanse |5| = r ~ ^/n, this takes 0{n) time. The final step involves 
inverting a x matrix for each i G V'. 

R is known that if r balls are placed in r bins uniformly at random, the maximum number 
of balls in any bin is lo'giogr 4“ ®(1)) with high probability (see [l9], for example). In this 
case qi = o(log r) for each point i , and hence the inversions take o(r log^ r) = o(n) time. 
Combining these facts, the algorithm runs in 0{Nlogn) time with high probability. □ 
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